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Abstract 

Background: The Fibroblast Growth Factor (FGF) pathway is driving various aspects of cellular responses in both 
normal and malignant cells. One interesting characteristic of this pathway is the biphasic nature of the cellular 
response to some FGF ligands like FGF2. Specifically, it has been shown that phenotypic behaviors controlled by 
FGF signaling, like migration and growth, reach maximal levels in response to intermediate concentrations, while 
high levels of FGF2 elicit weak responses. The mechanisms leading to the observed biphasic response remains 
unexplained. 

Results: A combination of experiments and computational modeling was used to understand the mechanism 
behind the observed biphasic signaling responses. FGF signaling involves a tertiary surface interaction that we 
captured with a computational model based on Ordinary Differential Equations (ODEs). It accounts for FGF2 binding 
to FGF receptors (FGFRs) and heparan sulfate glycosaminoglycans (HSGAGs), followed by receptor-phosphorylation, 
activation of the FRS2 adapter protein and the Ras-Raf signaling cascade. Quantitative protein assays were used 
to measure the dynamics of phosphorylated ERK (pERK) in response to a wide range of FGF2 ligand concentrations on 
a fine-grained time scale for the squamous cell lung cancer cell line H1703. We developed a novel approach 
combining Particle Swarm Optimization (PSO) and feature-based constraints in the objective function to calibrate 
the computational model to the experimental data. The model is validated using a series of extracellular and 
intracellular perturbation experiments. We demonstrate that in silico model predictions are in accordance with 
the observed in vitro results. 

Conclusions: Using a combined approach of computational modeling and experiments we found that 
competition between binding of the ligand FGF2 to HSGAG and FGF receptor leads to the biphasic response. At 
low to intermediate concentrations of FGF2 there are sufficient free FGF receptors available for the FGF2-HSGAG 
complex to enable the formation of the trimeric signaling unit. At high ligand concentrations the ligand binding 
sites of the receptor become saturated and the trimeric signaling unit cannot be formed. This insight into the 
pathway is an important consideration for the pharmacological inhibition of this pathway. 

Keywords: FGF signaling pathway, HSGAGs, Biphasic response, High throughput quantification, ODE-modeling, 
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Background 

Signaling pathways are arguably one of the most import- 
ant components driving how biological systems respond 
to environmental cues [1,2]. The ability of cells to per- 
ceive and respond to their microenvironment is the basis 
of development, tissue repair, and immunity as well as 
normal tissue homeostasis. Errors in cellular information 
processing contribute towards diseases such as cancer. 
By understanding the intricacies of cell signaling and 
processing, diseases may be treated more effectively. 

The FGF pathway plays a pivotal role in the stimulation 
of cell proliferation, cell migration and differentiation of a 
large number of cell types including muscle, neurons, cartil- 
age and bone cells [3,4]. FGF ligand - receptor signaling is 
regulated both by primary sequence differences between 
the 18 FGF ligands, the 7 main FGF receptors (FGFRlb, 
FGFRlc, FGFR2b, FGFR2c, FGFR3b, FGFR3c and FGFR4), 
by temporal and spatial expression patterns of FGFs, FGFRs 
and HSGAGs and FGF Binding Proteins. Tissue-specific 
alternative splicing in the second half of Ig-like domain 3 
(D3) of fibroblast growth factor receptors 1-3 generates 
epithelial FGFRlb-FGFR3b and mesenchymal FGFRlc- 
FGFR3c splice isoforms. This splicing event establishes a 
selectivity filter to restrict the ligand binding specificity of 
FGFRb and FGFRc isoforms to mesenchymally and epithe- 
lially derived fibroblast growth factors (FGFs), respectively 
[5]. FGF Binding Proteins (FGF-BP) have been described to 
function as a chaperone molecule that can mobilize FGF lo- 
cally and present it to the FGF receptor. The FGF-BPs have 
been described to enhance proliferation and signaling in 
NIH-3T3 cells [6]. 

HSGAG has been assigned multiple roles: it is known 
to serve as a co-receptor essential for signaling, as trans- 
port mediator to increase the local concentration of 
growth factors close to the cell surface or as a regulator 
to accelerate or attenuate signaling [7]. HSGAG are 
known to be essential for FGF signaling [8] and typical 
HSGAG levels on the cell surface are with 10 5 -10 6 mole- 
cules per cell [9] much higher than typical FGFR levels 
(< 50,000 molecules per cell Merrimack unpublished 
data). Thus, the benefits of understanding this pathway 
and the role of HSGAG in regulating FGF signaling are 
several fold: reveal greater insights into the fundamental 
principles of signaling pathway regulation by HSGAG in 
the case of FGF and more generally to understand the 
effect of inhibitors targeting the FGF pathway that are 
currently in development [10,11]. 

Disregulation of the FGF pathway components lead 
to various diseases, including multiple forms of malig- 
nant cancers [12]. FGFs are expressed ectopically in al- 
most 20% of different identified breast cancer cell lines 
[13]. It has also been shown that FGFs act as angio- 
genic growth factors that control capillary endothelial 
cell proliferation for vascular development [14]. A central 



paradigm of signaling pathways is that ligands bind to and 
activate cell-membrane bound receptors, which in turn 
leads to activation of intracellular cascades [15,16]. Typic- 
ally, binding of a monomeric ligand to a monomeric 
receptor follows Michaelis-Menten reaction kinetics. In- 
creasing the concentration of the ligand leads to an in- 
crease in receptor activation until ligand concentrations 
are high enough such that receptor activation is saturated. 
Intracellular receptor activation is followed by a cascade 
of enzymatic reactions that lead to the phosphorylation of 
effector molecules like ERK and AKT. Thus, increasing 
ligand concentration from low to intermediate levels 
increases the activation of ERK and AKT while at high lig- 
and concentrations, ERK and AKT are maximally acti- 
vated and therefore don't respond to further increases in 
ligand levels. Accordingly, one would expect that cells 
would respond to ligands in a saturable fashion as well 
(Figure 1A). This is indeed true for various signaling path- 
ways and physiological ligand concentrations, including 
ErbB and IGF1-R [17,18]. 

Interestingly, some cells respond to activation by FGF 
ligands atypically [3,4,19]. Instead of the typical satur- 
able response (Figure 1A), these cells respond in a bi- 
phasic manner (Figure IB). Specifically, cellular response 
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Figure 1 Schematic of sigmoidal (A) and biphasic (B) response. 
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increases from low to intermediate levels of FGF ligand 
but then decreases at high levels of the ligand. For in- 
stance, Williams et al demonstrated that FGF2-induced 
neurite outgrowth reaches a maximum at intermediate 
concentrations of FGF2 and that the outgrowth decreases 
at higher levels of FGF2 [3]. Garcia-Maya et al demon- 
strated that FGF-induced proliferation of NIH3T3 cells 
also reaches a maximum at intermediate concentrations 
of FGF2 [4,20]. Similar results have been shown for fibro- 
blasts and osteoblasts as well [19]. Notably, biphasic re- 
sponse to FGF ligands has most commonly been reported 
at the level of cellular phenotype, while the underlying 
molecular signaling events that lead to the biphasic re- 
sponse remain unexplored. 

Recently, Zhu et al. measured protein levels using 
Western Blots and indicated that FRS2- (FGFR adaptor 
protein) and ERK-phosphorylation is biphasic in re- 
sponse to FGFR activation by FGF2 [19]. Thus, they pro- 
vide hints that FGF signaling response does not follow 
the typical Michaelis-Menten reaction kinetics. How- 
ever, owing to the labor-intensive and low throughput 
nature of the experimental technique, protein levels 
were measured at a small number of time points and 
only for a small number of ligand concentrations. 
Moreover, they did not provide any mechanistic insight 
into how the FGFR pathway activation cascade drives 
this biphasic behavior. 

Efforts to model the FGFR pathway previously have 
primarily focused on the interactions taking place on 
the surface of the cell while completely ignoring the 
intracellular cascade [21-24]. Their efforts were di- 
rected towards understanding the effect of adding hep- 
arin, a soluble source of HSGAGs, on FGFR pathway 
activation. However, the time-course and dose-re- 
sponse of pERK to activation by ligand alone remains 
unexplored. In contrast, Yamada et al. built a combined 
mathematical model for the extracellular and intracel- 
lular components of the FGFR pathway but completely 
ignored the biphasic nature of response upon activation 
by the ligand [25]. 

Here, we investigate the mechanistic rationale for bi- 
phasic response to FGFR pathway activation utilizing a 
combination of high-density signaling data and Ordinary- 
Differential-Equation (ODE)-based mathematical model- 
ing. We report the dynamic pERK response to FGF2 
stimulation using a fine-grained time grid over 120 minutes. 
The model was calibrated to the experimental results using 
a published particle swarm optimization (PSO) algorithm 
[26] combined with a novel feature-based constrained 
optimization. Finally, we apply Markov-chain Monte Carlo 
sampling to explore nonidentifiability and uncertainty in 
the model calibration [27]. The computational model rep- 
resents the extracellular and intracellular components as 
well as the feedback regulation of the pathway. To validate 



the model, multiple simulations were conducted to predict 
signaling response to extracellular and the intracellular per- 
turbations of the pathway. The model predictions were val- 
idated by in vitro experiments. We demonstrate that 
without any fitting, the model explains each of these per- 
turbation experiments. We demonstrate that the complex 
protein interactions at the cell surface are necessary to ex- 
plain the observed biphasic dose-response while the nega- 
tive feedback loop from pERK to FRS2 controls the 
magnitude of the biphasic response. At low to intermediate 
concentrations of FGF2 there are sufficient free FGF recep- 
tors available for the FGF2-HSGAG complex to enable the 
formation of the trimeric signaling unit. At high ligand 
concentrations the ligand binding sites on FGFR become 
saturated and the trimeric signaling unit cannot be formed 
because binding of FGF2-HSGAG is weak, thereby leading 
to a decrease in pERK response. 

Results and discussion 

To uncover the underlying mechanism of biphasic FGF 
signaling response and to simplify the interpretation of 
the results, we use a representative non-small cell lung 
cancer cell line NCI-H1703 that was previously shown to 
primarily express FGFRlc and to induce Erkl/2 phos- 
phorylation upon stimulation with FGF2 [28]. The expres- 
sion of FGFRlc was confirmed by qPCR (Materials and 
methods Section 1). To calibrate the kinetic parameters in 
the model, dose-time matrices for ERK1/2 phosphoryl- 
ation in response to FGF2 were measured (Figure 2). A 
transient peak in ERK1/2 phosphorylation was observed 
at around 5 minutes of FGF2 exposure for all six doses. 
However, the rates of decay of ERK1/2 phosphorylation 
differed between doses, with sustained ERK1/2 phosphor- 
ylation observed even 120 minutes after ligand addition 
for intermediate doses (Figure 2B). This resulted in the bi- 
phasic dose response to FGF2 stimulation observed con- 
sistently after 10 minutes of exposure (Figure 2C). 

These detailed measurements increase the identifiabil- 
ity of the system, and allows the signaling cascade to be 
modeled mechanistically using a set of Ordinary Differ- 
ential Equations (ODEs) describing the biochemical re- 
actions. We constructed a mathematical description of 
the FGFR pathway that includes all the ligand-receptor 
interactions on the cell surface as well as the intracellu- 
lar MAP kinase signaling cascade. 

Description of the computational model 

The mathematical description consists of two major com- 
ponents - network topology and the corresponding 
network parameters. A schematic representation of the 
model structure is shown in Figure 3. The proposed struc- 
ture is based on previously published crystal structures 
with a 1:1:1 stoichiometry [29,30]. Briefly, the ligand binds 
to HSGAG to form a ligand-HSGAG complex, which 
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Figure 2 Mathematical model for pERK response to FGF2: the mathematical model explains the dynamics of pERK response to FGF2 
concentrations varying over three orders of magnitude. Experimental results are plotted as circles with standard deviations and model fits 
are plotted as solid lines. A). Feature-based constraints for fitting the mathematical model to pERK dose response data. B). Dynamics of pERK 
response to varying levels of FGF2-ligand (0.16 ng/ml to 500 ng/ml). pERK is measured at 1,2,3,4,5,8,10,20,30,60 and 120 min post ligand 
stimulation. C). pERK dose response curves at 5, 20 and 60 min post ligand stimulation. 



further binds to FGFR to form a trimeric complex [31]. 
This particular order of binding reactions was chosen for 
the formation of the trimeric complex since previous mea- 
surements have shown that FGF2-FGFR binding is much 
weaker than the binding of FGF2 to HSGAG. Moreover, 
they also showed that FGFR binding to to FGF:HSGAG 
complex is much stronger than FGFR binding to FGF2 
alone [32]. These results suggests that HSGAGs might 
have a regulatory function to control signaling by concen- 
trating FGF ligand on the cell surface, explaining how low 



concentrations of FGF are capable of activating signaling 
[7]. While other orders of binding reactions are possible, 
they are less likely to be important for the formation of 
trimer FGF:FGFR:HSGAG. Therefore, they have not been 
considered in this model to keep the size of the model 
and number of unidentifiable parameters small. The tri- 
meric complex dimerizes to form a 2:2:2 signaling unit 
that activates the intracellular signaling cascade. Specific- 
ally, the signaling unit binds to FRS2 and enzymatically 
phosphorylates it. In our model we assumed that FRS2 
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Figure 3 Schematic of the FGFR signaling pathway. Based on published surface plasmon resonance data, we assume that FGF binds first to 

HSGAG before it binds to the receptor. The Ras-Raf signaling cascade is reduced to a two-step cascade to increase model identifiability. The 

model also accounts for negative feedback from pERK to FRS2 and nuclear localization of pERK. 
k J 



directly binds to the signaling unit which is distinct from 
a model described previously [25]. pFRS2 subsequently ac- 
tivates the Ras-Raf signaling cascade. In the model, the 
intracellular signaling cascade is reduced to a two-step 
phosphorylation cascade. pFRS2 enzymatically phosphory- 
lates MEK and finally pMEK acts as an enzyme for phos- 
phorylation of ERK. The model also accounts for the 
negative feedback from pERK to FRS2 as described previ- 
ously: pERK binds to FRS2 and pFRS2 which ultimately 
leads to their degradation [33]. We also accounted for 
subcellular localization of pERK by considering its accu- 
mulation into the nucleus (detailed equations in Materials 
and methods Section 2). 

It should be recognized that the reactions included in 
this scheme are not necessarily direct protein-protein in- 
teractions. Some of these reactions represent a combin- 
ation of multiple interactions within the cascade. The 
underlying principle of such a reductionist approach is 
that for large non-identifiable systems like signaling path- 
ways, most components remain unmeasured. Therefore, a 
reduced mathematical model with appropriate features 
can better describe the experimental results while allowing 
for easy interpretation as compared to a model that ac- 
counts for each protein-protein interaction individually. 
The applicability of such reduced models is tested by 
verifying whether it can quantitatively fit experimental re- 
sults (Figure 2) with physically meaningful rate constants 
followed by validation of independent perturbation experi- 
ments without any fitting. 

Most of the kinetic rate constants described in the 
mathematical model remain unmeasured to date. Even 



for those rate constants that have been measured, for 
instance the ligand-receptor binding constants, a wide- 
range of values have been reported in the literature de- 
pending upon the experimental system and technique 
used [22,32,34]. Therefore, model calibration against 
the experimental data with physically meaningful pa- 
rameters is critical in order to gain a better understand- 
ing of the system. 

Calibration of the computational model using particle 
swarm optimization and feature-based constraints 

As is typical for signaling pathways, the mathematical 
model is parametrically non-identifiable with many more 
components than the number of measurements avail- 
able. Accordingly, there exists no unique solution for the 
parameter values and purely deterministic parameter es- 
timation techniques will fail to provide good estimates 
unless the initial parameter guess is in the vicinity of the 
global optimum itself. Therefore a global parameter esti- 
mation technique, Particle Swarm Optimization (PSO), 
was utilized to estimate the parameters of the model 
[35]. The approach allows for a more exhaustive and ef- 
ficient exploration of the objective function manifold to 
find good parameter fits. The version of PSO utilized in 
this report has been published previously and was par- 
ticularly useful for fitting parameters to signaling models 
[26,35]. The only inputs required for this approach are 
the physical constraints for the parameter values. 

Given the network structure, a critical aspect of esti- 
mating the parameters is the choice for the objective 
function to be minimized. We noted that minimization 
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of the most commonly used objective function, root 
mean square error (RMSE) between experimental data 
and simulation predictions was insufficient and the esti- 
mated parameters provide poor solutions. Therefore we 
developed a hybrid objective function that accounts for 
RMSE as well as differences between qualitative features 
of experimental and simulation data. Figure 2A summa- 
rizes the feature-based constraints used in the objective 
function. The objective function was penalized if the 
pERK level reached peak values in response to ligand 
stimulation too early or too late (detailed penalty func- 
tion expressions in Materials and methods Section 3). 
Similarly, the objective function was penalized if the 
pERK dose response at the relevant time points was not 
qualitatively biphasic. Minimizing the hybrid objective 
function using PSO led to physically meaningful param- 
eter estimates and was used to obtain the parameters 
used in this report. In contrast, using RMSE as the ob- 
jective function provided parameter estimates that had 
low objective function values but provided biophysically 
meaningless solutions. A detailed comparison between 
the two objective functions used and the corresponding 
parameter estimates obtained is described in Materials 
and methods Section 3. Based on these results, we 
propose a more general hypothesis that optimization of all 
signaling pathway models might benefit from utilizing a 
combination of qualitative constraints and RMSE values 
as compared to the simple RMSE-based objective func- 
tion. The applicability of this approach to other models re- 
mains to be tested and is the subject of future research. 

Using PSO and a combination of RMSE together with 
qualitative constraints as the objective function, multiple fits 
to the experimental data were estimated. Figure 2B shows 
the fit from one such typical data set in red and the gray 
lines represent fits of 6 additional parameter sets that are 
reported in Table 1. Additional file 1: Figure S1A shows that 
pERK response at all time-points is highly correlated (Pear- 
son correlation coefficient > 0.95) with the simulation results 
for multiple parameter sets. However, nonidentifiability of 
model parameters and uncertainty in parameter space in 
general can compromise the predictability of a mathematical 
model [36]. Therefore, parameter uncertainty was analyzed 
using the well-established Markov-chain Monte Carlo 
(MCMC) sampling approach. Recently, Hug et. al. showed 
that the approach was able to map uncertainty efficiently 
even for high-dimensional and non-linear models [27]. We 
tested the uncertainty of FGF model-parameters using a 
similar approach (detailed description in Materials and 
methods Section 6) and validated that the seven parameter 
sets identified by PSO are representative for the acceptable 
fits determined by the MCMC approach (Additional file 2: 
Figure S2). 

The model recapitulates pERK activation data at all 
different FGF2 concentrations. The model captures the 



experimental observations that at low FGF2 concentra- 
tions (below 4 ng/ml), the peak level of ERK activation 
increases with FGF2 concentration; however, at higher 
concentrations (above 4 ng/ml), the peak level of ERK 
activation remains constant (Figure 2B). Careful observa- 
tion reveals that at high FGF2 concentrations, the time 
to peak ERK activation decreases as FGF2 levels in- 
crease. Most importantly, the model also captures the 
fact that pERK de-activation varies with FGF2 concen- 
tration. Specifically, at intermediate levels of FGF2, 
pERK levels reach a peak after 5 minutes and then slowly 
decrease over the next hour. In contrast, at the high levels 
of FGF2, pERK levels peak before 5 minutes and decrease 
to low levels within 20-30 minutes. Therefore the dose re- 
sponse curve at time points later than 10 minutes shows a 
strong biphasic response (Figure 2C). Thus, the reduced 
model adequately captures all of the essential features of 
pERK dynamics in response to FGF2 stimulation. 

As described previously, even with the fine-grained 
pERK response measurements, the model is highly uniden- 
tifiable and therefore multiple parameter sets that might 
correspond to different response mechanisms explain the 
data equally well (Additional file 1: Figure S1A). Accord- 
ingly, before utilizing the aforementioned model and par- 
ameter sets shown in Figure 2 for further investigation of 
the underlying mechanisms, we tested and validated the 
model against multiple perturbation experiments. The val- 
idation will help rule out parameter sets that fit the training 
data-set well but do not describe the FGFR signaling path- 
way. Specifically, we perturbed the extracellular and intra- 
cellular components of the model and tested model 
validity by comparing the predictions with in vitro experi- 
mental measurements of pERK. The objective was to valid- 
ate the underlying mechanisms predicted by the model 
rather than exact quantitative numbers. Therefore, the 
model was not trained on any these new experimental re- 
sults and all model predictions were qualitatively compared 
to the experimental results. 

Addition of soluble heparin decreases the pERK response 
at low and intermediate FGF2 levels but increases the 
pERK response at high levels of FGF2 

One of the essential components of the model is the 
binding of FGF2 to HSGAGs. It is known that FGF2 also 
binds to soluble heparin and this binding in solution can 
compete with the binding of FGF2 to HSGAG. Soluble 
heparin can also rescue signaling behavior in cells that 
have been stripped of cell surface HSGAGs. Therefore, 
to validate the ligand binding and signaling complex for- 
mation module of the model, we tested how the pERK 
response changes with addition of soluble heparin to 
media. The effect of heparin on the FGF pathway has 
been explored previously at the level of surface interac- 
tions and ligand-receptor complex formation [21,24]. 
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Table 1 Parameter values 



Name 


Unit 


Setl 


Set2 


Set3 


Set4 


Set5 


Set6 


Set7 


Fixed parameters 
















N A 


molecules/mole 


6.02E+23 


6.02E+23 


6.02E+23 


6.02E+23 


6.02E+23 


6.02E+23 


6.02E+23 


Vshell 


liter 


4.85E-15 


4.85E-15 


4.85E-15 


4.85E-15 


4.85 E- 15 


4.85 E- 15 


4.85E-15 


Ncell 


cells 


4.00E+04 


4.00E+04 


4.00E+04 


4.00E+04 


4.00E+04 


4.00E+04 


4.00E+04 


Vfluid 


liter 


1 .24E-04 


1 .24E-04 


1 .24E-04 


1 .24E-04 


1 .24E-04 


1 .24E-04 


1 .24E-04 


FGFR 


molecules/cell 


2.00E+04 


2.00E+04 


2.00E+04 


2.00E+04 


2.00E+04 


2.00E+04 


2.00E+04 


Fitted parameters 
















HSGAG 


molecules/cell 


1 .OOE+05 


1 .OOE+05 


1 .27E+05 


1 .OOE+05 


1 .OOE+05 


1 .OOE+05 


1 .OOE+05 


FRS2 


molecules/cell 


1 .OOE+04 


5.27E+05 


4.37E+04 


1 .OOE+04 


1 .OOE+04 


1.16E+04 


1 .02E+04 


MEK 


molecules/cell 


1 .00E+06 


7.24E+05 


3.08E+05 


3.20E+05 


3.73E+05 


1.55E+05 


3.66E+05 


ERK 


molecules/cell 


1 .30E+06 


1 .52E+06 


1 .34E+06 


1 .88E+06 


4.93E+06 


1 .09E+06 


3.61 E+06 


kfo 


1/((molecule/cell)*s) 


2.67E-08 


1.10E-08 


8.61 E-09 


1 .97E-08 


1 .46E-08 


1 .52E-08 


2.60E-08 


Ko 


1/s 


1 .94E-03 


8.00E-04 


6.29E-04 


1 .44E-03 


1 .07E-03 


1.11 E-03 


1.95 E-03 


kfla 


1/((molecule/cell)*s) 


7.89E-10 


1.13E-09 


3.14E-09 


9.28E-10 


1 .32E-09 


2.04E-09 


2.37E-09 


k r ia 


1/s 


9.17E-05 


1.31 E-04 


3.63E-04 


1 .07E-04 


1 .52E-04 


2.36E-04 


2.75 E-04 


kf5a 


1/((molecule/cell)*s) 


8.62E-09 


1 .72E-08 


2.47E-08 


5.65E-09 


1.19E-08 


6.57E-08 


3.10E-06 


k r 5a 


1/s 


2.60E-07 


1 .27E-03 


1 .90E-04 


3.03E-08 


2.84E-06 


3.24E-07 


2.14E-01 


kfdim 


1/((molecule/cell)*s) 


8.48E-02 


1.50E-01 


2.05 E-03 


1.30E-01 


2.01 E-02 


5.26E-02 


3.87E-06 


krdim 


1/s 


9.94E+00 


1 .OOE-05 


2.90E-04 


1 .92E-04 


1.09E-01 


5.83E-04 


1.35 E-02 


kfph 


1/s 


3.30E+00 


4.92E-01 


1 .58E-02 


5.36E+00 


2.49E-01 


5.02E-02 


2.08E-03 


kfintl 


1/s 


1 .38E-04 


2.14E-03 


1 .00E-06 


9.53E-05 


9.90E-06 


4.38E-04 


2.45 E-05 


krintl 


1/s 


4.70E-07 


4.26E-04 


2.52E-03 


2.31 E-07 


2.42E-03 


9.73E-08 


9.03E-04 


kfis 


1/((molecule/cell)*s) 


2.41 E-04 


4.61 E-04 


5.86E-04 


2.09E-04 


2.14E-04 


4.43 E-04 


1.73 E-04 


k r i5 


1/s 


4.21E-03 


4.87E-06 


3.40E-03 


1 .66E-04 


2.11 E-02 


8.33E-05 


1 .20E-05 


km 


1/s 


4.72E-01 


5.47E-01 


5.88E+00 


3.95E+00 


1.00E+01 


5.50E-01 


9.03E+00 


kf 3 5 


1/((molecule/cell)*s) 


3.58E-04 


5.43E-04 


1 .48E-04 


1 .47E-04 


9.23E-05 


1 .40E-04 


5.80E-05 


k r 35 


1/s 


2.34E-01 


1.13E-01 


5.24E-05 


2.74E-02 


2.00E-03 


4.1 5 E-06 


1 .70E-05 


kf36 


1/s 


3.04E-01 


2.34E-01 


5.53E-02 


1.33E-02 


2.57E-02 


4.08E-02 


3.78E-02 


kf 3 9 


1/((molecule/cell)*s) 


2.61 E-05 


4.14E-08 


3.28E-07 


8.1 5 E-06 


1 .00E-03 


8.73E-05 


5.78E-06 


k r39 


1/s 


1 .66E-06 


4.91 E-06 


6.60E-05 


3.94E-05 


2.99E-01 


3.55E-04 


1 .04E-04 


kf40 


1/s 


9.95E-02 


2.62E-01 


1.11E-01 


1 .74E+00 


1.51 E+00 


1.71 E-01 


1 .04E+00 


kfdpi 


1/s 


4.71 E+00 


6.65 E+00 


4.29E+00 


9.96E+00 


1.00E+01 


5.43E+00 


5.45 E+00 


kfd P 2 


1/s 


1 .03E+00 


8.93E-03 


1.83E-02 


3.10E-02 


8.29E+00 


1 .78E+00 


6.56E-02 


kfd P 3 


1/s 


1 .28E+00 


6.37E-03 


5.71 E-03 


4.66E-03 


4.64E-03 


5.07E-03 


1.51 E-02 


kf43 


1/((molecule/cell)*s) 


1 .06E-04 


9.66E-06 


3.08E-06 


3.13E-07 


1.53E-07 


8.09E-07 


3.92E-07 


k r43 


1/s 


3.34E-05 


7.34E-01 


5.43E-01 


1.71 E-04 


1.71 E-04 


2.21 E-02 


8.90E-05 


kf44 


1/s 


2.14E-04 


4.28E-02 


2.10E-01 


4.39E-04 


3.91 E-04 


3.05E+00 


3.05E-04 


kf47 


1/s 


3.36E+00 


2.32E-04 


9.67E-05 


2.10E-05 


2.01 E-05 


1 .52E-04 


8.06E-02 


k r47 


1/s 


1 .47E-02 


8.25 E-05 


6.01 E-07 


1.15E-07 


9.96E-02 


5.54E-06 


4.00E-02 



However, those results and interpretations have never been 
extended to include the effect on the intracellular signaling 
cascade. Using our model, we explore the effect of heparin 
on pERK dynamics at the same fine time-grid as discussed 
in the previous section (Materials and methods Section 4). 



Model simulation of heparin addition to extracellular 
medium provided some non-intuitive insights (Figure 4A, 
B solid curves). At low to intermediate levels of FGF2, 
addition of heparin was predicted to decrease the level of 
ERK activation at all the time points. This is in line with 
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Figure 4 Model validation: extracellular pathway perturbation. A). pERK time-response curves up to 2 hrs upon addition of FGF2 with (red) 
and without (black) soluble heparin. Experimental results are plotted as circles with standard deviations and model fits are plotted as solid 
curves. Note that the red curves are model predictions and not fitted to pERK response to FGF2 with heparin. B). Model prediction for pERK 
dose-response curve at 10 min upon addition of varying amounts of heparin into the extracellular medium. C). Experimental validation of pERK 
dose-response upon heparin addition to the medium. 



the role of heparin as a ligand trap. However, surprisingly, 
at high levels of FGF2, addition of heparin was predicted to 
increase the level of pERK response. Model prediction for 
change in time-course and dose-response of pERK upon 
addition of heparin were validated using in vitro experi- 
ments (Figure 4A, C symbols and dotted curves). 

Without any fitting, the model accurately captures 
pERK response to FGF2 in the absence/presence of hep- 
arin (Figure 4). Note that the same qualitative predic- 
tions were made by all the different parameter sets 



(Additional file 1: Figure SIB). Therefore, the level of 
confidence in the explanation provided by the simula- 
tions is high: The model indicates that at low levels of 
FGF2, heparin acts as a ligand trap and reduces the level 
of FGF2 binding to HSGAG and subsequent formation 
of signaling complexes while at high levels of FGF2, 
there is enough excess FGF2 present in the extracellular 
medium to overcome the ligand trap and to form signal- 
ing complexes. Additionally, a small number of heparin- 
FGF2 complexes bind to FGFR and initiate signaling 



Kanodia et al. Cell Communication and Signaling 2014, 12:34 
http://www.biosignaling.eom/content/1 2/1 /34 



Page 9 of 1 8 



complex formation, leading to an increase in pERK re- 
sponse. Thus, without fitting any of the parameters for 
soluble heparin, the model captures pERK response to 
extracellular perturbations. It is noteworthy here that 
some experimental data points do not match perfectly 
with simulation results. We hypothesize that the differ- 
ences could be due to non-pathway specific interactions 
or other higher order interactions that are indeed a part 
of the FGFR pathway but have a small influence on the 
overall response and are not captured by the current 
model. Explanation of these differences would require 
measurement of more proteins within the cascade com- 
bined with a more detailed model and will be the subject 
of future work. 

Delayed inhibition of MEK post FGF2 stimulation 
amplifies the biphasic pERK response 

To validate the intracellular ERK activation module of 
the pathway, MEK was inhibited using a small molecule 
inhibitor U0126 (henceforth referred to as 'MEKi') with 
two different schedules. In the first experiment, MEKi was 
added to the media at the same time as the ligand itself 
(Figure 5A). Five different concentrations of MEKi were 
used in combination with seven different concentrations 
of FGF2. For each combination of MEKi-FGF2, the pERK 
response was measured in cells 10 minutes after FGF2 
stimulation. As predicted by the computational model 
only the 0.37 \iM and 1.11 \iM concentrations of MEKi 
lead to a significant pERK inhibition across the entire 



range of FGF2 concentrations. It is noteworthy that the bi- 
phasic shape of the FGF2 dose-response curve remained 
the same for different MEKi levels. pERK levels reach a 
peak at 4-20 ng/ml of FGF2 and then decrease at higher 
FGF2 concentrations. These results were in line with 
the predictions made by the model. The ODEs used to 
describe MEKi interactions are outlined in Materials 
and methods Section 5 in the Appendix. 

It is not surprising that co-incubation of MEKi with 
FGF2 was expected to show uniform inhibition of pERK. 
Therefore, another experiment was conducted where the 
MEK inhibitor was added 5 minutes post ligand stimula- 
tion and pERK response was measured 5 minutes post 
MEKi addition (Figure 5B). In this experiment, the cells 
respond to ligand stimulation for 5 minutes under 
MEKi-free conditions and then for 5 minutes in the 
presence of MEKi. Thus, this perturbation experiment is 
a more stringent test of the model since it is required to 
accurately predict the combination of control as well as 
perturbation conditions within the same simulation. 
Similar to the co-inhibition experiment in Figure 5A 
the model predicted reduced pERK levels - however to 
a much lesser extent than in the co-inhibtion experi- 
ment. Specifically, the model predicted that the peak 
ERK response at intermediate FGF concentrations is sig- 
nificantly inhibited while the signal at the low and high 
FGF concentrations are inhibited to a lesser extent. This 
prediction was validated as showed in Figure 5B in the 
delayed addition experiment. The fact that the model 




0.1 1 10 100 

FGF2 [ng/ml] 



B 



FGF 
l_ 



5 min 



CD 

^ 1.0 
co 

E 0.8 
g 0.6 



M ^ 0.4 
5 min . DC 

LU 0.2 
Q. 

0 



MEKi 

J_ 



I 

pERK 




0.1 1 10 100 

FGF2 [ng/ml] 



1.2 
1.0 
0.8 
0.6 
0.4 
0.2 



V-tT 



0.1 1 10 100 

FGF2 [ng/ml] 



1.2 
1.0 
0.8 
0.6 
0.4 



0.2ft 

0 



J x V 

Ay * 



0.1 1 10 100 

FGF2 [ng/ml] 



- no inhibitor - MEKi = 0.75*N(MEK) ■■■ no inhibitor ■■■ 0.37 \xM MEKi 

- MEKi = 0.25*N(MEK) - MEKi = 1*N(MEK) ... 0.12 pM MEKi ■■■ 1.11 nM MEKi 

Figure 5 Model validation: intracellular pathway perturbation. A). pERK dose response curve at 10 min post ligand stimulation for 
simultaneous addition of FGF2 ligand and MEKi model prediction compared to experimental validation. B). pERK dose response curve at 10 min 
post ligand stimulation for MEKi addition 5 min post FGF2 ligand stimulation model prediction compared to experimental validation. 
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captures the partial inhibition of pERK in the delayed 
inhibition experiment indicates that the negative feed- 
back loop from pERK to FRS2 is accurately captured by 
this parameter set. It should be noted here that there 
are differences between pERK response under control 
conditions (i.e. no MEKi) for the two experimental setups. 
This difference can be attributed to the super-sensitivity 
of pERK to fluid movement and addition rather than 
biological response to growth factors or inhibitors. There- 
fore, it is further emphasized in this case to focus on 
the qualitative prediction of model compared with ex- 
perimental results and interpret the changes in pERK 
response with respect to control conditions rather than 
absolute values. 

It is also noted that not all the parameter sets pre- 
dicted the effect of MEKi accurately. One of the sets 
gave almost no response to MEKi while two others in- 
correctly predicted pERK response to stimulation by 
FGF2. As expected previously, matching model predic- 
tion with the more stringent experimental test of fitting 
pERK response to FGF2 and MEKi scheduled with a 
time lag helps identify physiologically relevant parameter 
sets for the FGFR pathway. Finally, only four out of the 
seven parameter sets that predicted both perturbation 
experiments correctly were used for subsequent model 
interpretation and analysis. Overall, one can be fairly 
confident that the model significantly captures the be- 
havior of intracellular ERK activation by FGF pathway 
and that the observed biphasic pERK dose response can 
be explained by the complex interplay between FGF re- 
ceptor, FGF2 ligand and HSGAG. 

Conclusion 

Biphasic phenotypic responses to FGF signaling have 
been observed in multiple cell types in varying contexts. 
However, an investigation of the mechanism behind bi- 
phasic response at the level of signaling components 
remains unexplored. Here, we utilize high-throughput 
measurements of ERK activation for measuring the dy- 
namic signaling response to a wide range of FGF2 con- 
centrations. Such detailed quantification allows us to 
tease apart subtle changes in ERK activation dynamics 
and thus provide a more complete picture of the system. 
The quantification also serves as an excellent training 
dataset for a mechanistic ODE-based model of the FGFR 
pathway. Given the lack of identifiability of the model, 
the global parameter estimation method PSO was uti- 
lized with a custom-built objective function to estimate 
the parameters. Finally, we validated that the identified 
sets of parameters are representative using an MCMC 
approach. We demonstrate that the model successfully 
captures the subtle changes in the timing of pERK peak 
levels as well as the change in pERK decay rates. To 
validate the assumptions made to build the pathway 



network structure and the estimated parameter values, 
the model was tested against multiple perturbation ex- 
periments. Without any fitting, the model accurately 
predicted changes in pERK response to both extracellu- 
lar and intracellular perturbation conditions. Thus, we 
assert that the reduced model captures the essential fea- 
tures of the system appropriately and can be utilized for 
further investigation. It is noteworthy that the objective 
of our approach was to build a mathematical model 
based on commonly understood aspects of FGFR path- 
way biology and uncover a plausible mechanism for 
biphasic pERK response. Therefore, we did not exhaust- 
ively explore the space of model structures for alterna- 
tive theoretical models that can also explain the data. 
Such an exploration can be instrumental for uncovering 
completely unknown or less well-understood biology 
and will be the subject of future research. 

One of the primary goals of the model is to help un- 
cover the underlying mechanism of biphasic signaling 
response to activation by FGF2 ligand. In our model we 
assume a 2:2:2 stoichiometry with FGF binding first to 
HSGAG and subsequent binding of the FGF:HSGAG 
complex to FGFR before dimerization occurs. A close 
look at the model indicates that the competition be- 
tween binding of FGF2 ligand to HSGAG and FGFR 
leads to the observed biphasic response. At low to inter- 
mediate concentrations of FGF2, despite the binding of 
ligand to both HSGAG and FGFR, there are enough free 
FGFRs on the cell surface for the FGF2-HSGAG com- 
plex to bind and initiate a trimeric signaling-unit forma- 
tion. However, at high levels of FGF2, ligand binding 
sites become saturated; specifically, a large fraction of 
the FGFR molecules are bound to FGF2 and trimeric 
signaling units cannot form, because binding of FGF2- 
HSGAG is weak, thereby leading to a decrease in pERK 
response (Figure 6A). Recently, Brown et al published 
data that favors a pre-assembly model where FGF binds 
first to HSGAG in a transdimeric manner before it binds 
to two receptors forming a 2:2:1 complex [37]. This 2:2:1 
structure aggravates the competition for free FGFRs to 
form a signaling unit compared to the 2:2:2 structure 
implemented in our model. This strengthens the model- 
based prediction that lack of free FGFRs at high ligand 
concentration can explain the observed biphasic behav- 
ior. Although the presence of biphasic response is driven 
by surface interactions, the strength of biphasic-ness is 
also regulated by the intracellular cascade. For instance, 
in the absence of pERK-pFRS2 feedback loop, pERK 
levels increase according to our model simulations in re- 
sponse to FGFR activation and remain maximal until the 
receptor is internalized and degraded. In this case, the 
time point at which pERK reaches peak level changes for 
each FGF2 concentration but the sustained biphasic re- 
sponse after 5min can not be observed. Therefore, the 
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Figure 6 Mechanism of biphasic response. A). Schematic of the model explaining why high levels of FGF2-ligand lead to a lower pERK signal 
as compared to intermediate levels of FGF2. B). pERK dose-response in a cell line with higher FGFR and HSGAG molecules (MDAMB-134) 
comparing model prediction to experimental validation. 



nature and the level of biphasic response are regulated 
by a combination of receptor competition on the surface 
as well as signaling feedback inside the cells. 

The computational model was explored in silico and 
then tested experimentally by using another cell-line. 
According to the model, a system with a large number 
of FGF receptors is predicted to show a sigmoidal pERK 
response (Figure 1A). At high FGF receptor levels, there 
is no competition between FGF2 binding to FGFRs and 
HSGAG and thus, as expected, pERK response saturates 
at high FGF2 levels rather than showing a biphasic 
response. This model interpretation was verified experi- 
mentally by using another cancer-cell line MDA-MB- 
134 (Figure 6B). This cell line has about 8x more FGFR1 
mRNA and ~2x more HSGAG and syndecan-1, a cell- 
surface proteoglycan containing HSGAGs [38]. Under 
high FGFR and high HSGAG conditions, the model 
predicts that the pERK response looks sigmoidal up to 
500 ng/ml of FGF2. The model prediction was verified 
by measuring pERK in MDA-MB-134 cells 10 min post 
FGF2 stimulation. These results highlight the import- 
ance of HSGAG and soluble heparin in regulating FGFR 
pathway activation. They indicate that HSGAG is not 
just a passive scaffold that facilitates binding of FGF 



ligands to FGFRs but actively participates in the local 
regulation of the dynamics of intracellular signal activa- 
tion [5]. It is noteworthy at this point that currently 
there are no good methods for accurately quantifying 
levels of HSGAGs. HSGAGs are not homogeneous, but 
rather a mixture of HSGAG species of different lengths 
and sulfation patterns. The sugar sequence and level 
of sulfation can impact the level of binding and its 
role in the signaling complex [39]. Therefore, in the 
absence of such measurements, the model provides an 
effective alternative approach for investigating the im- 
portance of HSGAGs in modulating signaling response 
to FGF ligands. 

The FGFR pathway signaling model and estimated par- 
ameter fit presented in this report is specific for the FGF2 
ligand interacting with the FGFR1 receptor. However, the 
quantitative approach, including measurements, model 
topology and estimation approaches are generalizable. 
HSGAGs affect the kinetics of ligand-binding and thus 
signaling of several other growth factors including HB- 
EGF, HGF, PDGF, MBPs, TGF-p or wingless (a member of 
the Wnt family) [40-45], In the case of HB-EGF binding 
to the EGF receptor, binding is enhanced at relatively 
low levels of heparin but higher heparin concentrations 
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lead to inhibition of receptor, as is the case in FGF sys- 
tem [42]. In contrast, despite strong in vivo evidence 
for a stimulatory role for endogenous HSGAGs in 
bone-morphogenetic-protein (BMP) signaling, the re- 
sults of in vitro studies have been inconsistent [46,47]. 
Since FGF and BMPs differ in the way they form signal- 
ing competent dimers, the biphasic response that can 
be observed with FGF may not be observable with 
BMPs or other growth factors [48]. Therefore, the role 
of HSGAGs in different signaling pathways will need to 
be evaluated individually. 

In summary, this work increases the quantitative un- 
derstanding of the FGF signaling pathway and its intri- 
cate regulation on the cell surface as well as downstream 
of the receptor, which is important to understand cellu- 
lar decision making. 

Materials and methods 

Section 1 . Experimental details of quantitative pERK 
measurements & perturbation experiments 
Cell culture 

NCI-H1703 cells were obtained from American Type Cul- 
ture Collection (ATCC CRL-5889). Cells were maintained 
in RPMI 1640 medium (Lonza) supplemented with 10% 
fetal calf serum (Tissue Culture Biologicals, Lot # 107197), 
2 mM L-Glutamine (Gibco), penicillin (100 U/mL) and 
streptomycin (100 (ig/mL) (Gibco), and grown in a hu- 
midified atmosphere of 5% C02 and 95% air at 37°C 

In vitro signaling studies 

For in vitro signaling studies, cells were seeded in 96-well 
tissue culture plates, allowed to attach overnight, then 
switched to serum-free media supplemented with 0.5% 
bovine serum albumin ( Sigma- Aldrich) for 16 to 20 hours. 
For dose-time matrices, serum-starved cells were stimu- 
lated with serial dilutions of FGF-2 (R&D Systems, 4114- 
TC-01 M) starting at 500 ng/mL for 1, 2, 3, 4, 5, 7, 10, 15, 
30, 60, 120 minutes. After stimulation, cells were washed 
with ice cold phosphate-buffered saline (PBS), then lysed 
in cold lx lysis buffer from the AlphaScreen SureFire p- 
Erkl/2 Assay Kit (Perkin-Elmer TGRES10K). 

For perturbation with exogenous heparin addition, cells 
were seeded as described above, then switched to serum- 
free media supplemented with 0.5% bovine serum albumin 
with or without 500 ug/mL heparin sodium salt (Sigma- 
Aldrich H3149-250KU). Serum-starved cells were stimu- 
lated with serial dilutions of FGF-2 in the starvation media 
and lysed as before. For perturbation with MEK inhibition, 
cells were seeded and starved as described above. Serum- 
starved cells were treated with one of two protocols: 1) 
simultaneous addition of a dose matrix of serially diluted 
FGF-2 starting at 500 ng/mL and U0126 starting at 10 uM 
(EMD Chemicals 662005-1MG) for 10 minutes followed 
by PBS wash and lysis, 2) addition of serially diluted 



FGF-2 followed by addition of serially diluted U0126 
after 5 minutes of incubation and finally PBS wash and 
lysis 10 minutes after ligand addition. 

Measurement of phospho-Erkl/2 levels on all samples 
were performed with the AlphaScreen SureFire pErkl/2 
assay kit (Perkin-Elmer TGRES10K), according to the 
manufacturer protocol in 384-well plate format (Alpha 
Plate, PerkinElmer 6008350) and read using the EnVision 
2013 Multilabl Plate Reader (Perkin-Elmer). 

FGF Receptor Expression by qPCR 

RNA was isolated from 2xlO A 6 cells using RNeasy kit (Qia- 
gen) following manufacturer protocols. Genomic DNA was 
eliminated using deoxyribonuclease (DNase) treatment 
using DNase I (Roche). 6 g RNA was reverse transcribed to 
cDNA in a 60 L reaction using a High Capacity cDNA Re- 
verse Transcription kit (Life Technologies Cat#4368814), 
and samples were stored at -80°C until qPCR. 

FGFR isoform expression levels were measured in 
duplicate in a 384-well reaction plate on an Applied Biosys- 
tems Viaa7 instrument (Life Technologies) with SYBR 
Green chemistry. 50 ng cDNA per well, and primer con- 
centrations of 150 nM were used. Primer pairs were de- 
signed to detect receptor Illb and IIIc isoforms selectively 
based on Fon Tacer et al. [49]. (see Table 2 below) For GusB 
and GAPDH, primers were purchased from Integrated 
DNA Technologies (Hs.PT.5 1.2648420, Hs.PT.51.1940505 
respectively). FGFR primer pairs were validated using plas- 
mids from Origene for specificity and efficiency. 

qPCR data was analyzed by Viaa 7 Software vl.2.2 . Base- 
line values were set automatically, and threshold values 
were kept constant. Samples with Ct values of >35 were 
considered below the limit of detection. Expression levels 
were normalized to the average of the housekeeping genes 
using the delta Ct method, and normalized expression is 
calculated as 2 A -deltaCt. These results are plotted below. 

Results - The relative expression of each of the FGF re- 
ceptor isoforms was measured in NCI-H1703 cells as part 
of a wider screen. Based on the expression levels, FGFRlc 
contributed more than 95% of the total while FGFRlb, 

Table 2 Primer sequences for FGFRs 



Primer name Sequence (5'-3') 



FGFR1 forward 


CAACCTGCCTOTGTCCAGATC 


FGFR1 Illb reverse 


CTCCGCATCCGAGCTATOA 


FGFR1 IIIc reverse 


ATCTCmGTCGGTGGTATOACTC 


FGFR2 forward 


GGGCTGCCCTACCTCAAG 


FGFR2lllb reverse 


GCCAGCAOTCTGCATOGA 


FGFR2lllc reverse 


ATCTCmGTCCGTGGTGT 


FGFR3 forward 


acggcacaccctacg™ 


FGFR3lllb reverse 


ACGTCGGCCTCCACACTCT 


FGFR3lllc reverse 


CTCOTGTCGGTGGTGTOGC 
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FGFR2b, FGFR2c, FGFR3b and FGFR3c combined con- 
tribute less than 5% of the total Consistent with previous 
reports, NCI-H1703 cells predominantly express FGFRlc 

Section 2. Equations and parameters (table of seven sets) 
of the mathematical model 

Additional file 3: Figure S3 summarizes the reaction net- 
work of the model in graphical form. The ordinary dif- 
ferential equations and parameters used are as described 
below (Table 3). 

Set of ordinary differential equations 

-{kfla * S 2 * Si-k r i a * S 9 )-(k f5a * Si * S 8 -k r5a * S W ) 
" (kfo * S 2 * S 3 -kr0 * S S ) ~ (kfia * 5 2 * S X -k r i * S 9 ) 
-(kfo *5 2 *S 3 -krQ *5 8 ) 



= ~{kf35 * Su * 5 5 -/< r35 * 5 2 l ) + kf dp 2 * « 



ds\ 
dt 

ds 2 
dt 



ds 3 
dt 



ds^ 
dt 



-(kf 15 * 5 4 * s n -k r is * 5 20 ) 

+kfdpl * 5i4-(/C/43 * «19 * H-k r ^ * 5 i6 ) 



Table 3 List of species 


Name 


Species 


Si 


FGFR 


s 2 


FGF 


S3 


HSGAG 


s 4 


FRS2 


s 5 


MEK 


s 6 


ERK 


s 7 


pERK nuc | ear 


s 8 


FGF:HSGAG 


s 9 


FGF:FGFR 


SlO 


FGF:HSGAG:FGFR 


S11 


(FGF:HSGAG:FGFR) 2 


Sl2 


p(FGF:HSGAG:FGFR) 2 


S13 


i(FGF:HSGAG:FGFR) 2 


S14 


pFRS2 


Sl5 


pMEK 


Sl6 


pERK:FRS2 


S17 


FRS2 u bj qu j t j natec | 


Sis 


pERK:pFRS2 


Sl 9 


pERK 


s 20 


p(F:H:R) 2 :FRS2 


S 2 1 


pFRS2:MEK 


S 2 2 


pMEK:ERK 



ds { 



<fdp2 * S15 

= -(kf 39 * 515 * s 6 -k r39 * S22) + kf dp3 * 5i9 
^ = (/C/47 * 5i9"/<r47 * 5 7 ) 

ds% 

(kfO * 52 * 5 3 -/c r0 * 5 8 ) - (kfsa * Si * 5 8 -/r r5a * s 10 ) 

-(*/l« *5 2 *Si-k r ia *5 9 ) 

(*/5« * 5l * S 8 -k r5a * 5io)-2 
* (/^m * 5i 0 * SiQ-krdim * 5n) 

= (*/Hww * 5i 0 * SiQ-krdim * $\\)-kfph * 5n 

-(*/2*fl * 5i2"/Wl * 5i 3 ) + A^fc * 5n 
+kf 19 * 520-(*/l5 * 512 * 5 4 -/<rl5 * 5 2 o) 

= -{kfintl * 5i2"/Wl * 5i 3 ) 



ds 9 

dt 

ds\Q 
dt 

ds 
d\ 

ds\2 
dt 

ds 



ds 



= k f i 9 * s 20 -kfdpi * 5i 4 -(/r /3 5 * 514 * s 5 -k r35 * s 2 i) 

+ kf 36 * 52l-(/C/43 * 519 * 5i 4 -/Cr43 * 5l 8 ) 



ds 



= kf 36 * S 2 \-kfdp2 * Si$-(kf 39 * 5i 5 * S 6 -I< r39 * 5 22 ) 
+ /C/40 * 5 2 2 



= (AT/43 * 519 * 5 4 -/<r43 * 5i 6 )-At/44 * $16 

= A:/44 * 5i6 + &/44 * 5i8 
^ = * 519 * 5i4-/<r43 * Sis) -kf44 * $18 



dSy] 

dt 
dsis 



ds 



= kfw * S 2 2-kfd p3 * 5i9-(/C/ 43 * 5 i9 * 5 4 -/c r43 * 5 i6 ) 
+/C/44 * 5i6-(A:/43 * 5l9 * 5i4-/<r43 * 5l 8 ) 
+/C/44 * 5i 8 -(/C/47 * 5l9-^r47 * 5 7 ) 



ds' 



-^jT = (kf 15 * 5 4 * Si2-k r lS * 5 2 o)-*/l9 * 5 2 0 
= (kf 35 * 514 * S 5 -k r35 * S 2 l)~kf36 * $21 

= (*/39 * 515 * S 6 -J< r39 * 5 22 )-/<:/40 * 5 2 2 



dS2l 

dt 

ds 2 2 
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Section 3. Comparing different forms of objective 
function for optimization 

One of the most commonly used forms of objective func- 
tion for minimization using local and global optimization 
methods is the root mean square error (RMSE). The 
error term serves as a scale-dependent estimation of the 
deviation between observed and predicted values. This 
functional form is the exact equation that needs to be 
minimized for fitting a linear relationship between input 
and output. However, for highly non-linear models, such 
as the ODEs described above, the RMSE-based objective 
function manifold in the multi-parameter space is rife with 
a large number of local minima and therefore not suitable 
for parameter estimation. An alternative to RMSE-based 
objective functions is to construct hybrid objective func- 
tions that combine RMSE error with feature-based con- 
straints incorporated as cost functions. The details of the 
hybrid objective function used for fitting the pERK data 
(Figure 2B) to FGFR model (Materials and methods, 
Section 2) are as follows. 

Based on the physics of biological systems, the follow- 
ing constraints were incorporated into the objective 
function. 

1. Total number of phosphorylated molecules of all 
intracellular molecules should exceed 100 at peak 
activation (This constraint needs to be satisfied so 
that the mean field approximation required for 
building ODE-based models is valid). If this was not 
fulfilled the following penalty was added 



^=nn 



100 



lll±m(*,;) 

where m is equal to the number of FGF concentrations 
used, n is equal to the number of intracellular molecules 
considered (ERK, MEK, FRS2 and FGFR) and m(i,j) is 
equal to the peak level of phosphorylated molecules of 
species j for FGF concentration i. 

2. pERK should be a sizable fraction of total-ERK 
present in the system (This constraint needs to be 
satisfied to ensure that the parameters don't fall 
into a regime where only a small fraction of pERK 
molecules get phosphorylated and control the 
system). 



of totERK that must be phosphorylated and max(pERK) 
is equal to the maximal level of pERK molecules across 
all doses of FGF treatment. 

3. The number of phosphorylated molecules down the 
cascade should not decrease by more than 25% 
within a single step (This constraint is expected to 
be true since kinases act as catalyst for 
phosphorylation of the molecule down the cascade 
and therefore it is expected that the input signal 
amplifies down the cascade). 



0.75 * max (kinase(ij)) 
max (substrate(ij)) 



where m is equal to the number of FGF concentra- 
tions, n+1 is equal to the number of kinases in the 
system (FRS2, MEK and ERK) and max(kinase(i,j)) or 
max(substrate(i,j)) is maximum level of upstream or 
downstream kinase respectively for FGF dose i and 
kinase 

Based on the experimental data, the following con- 
straints were incorporated into the objective function 

4. The pERK time response curve should reach a peak 
between 4-10 min after ligand stimulation 



if t(peak) > 10 min : ^4 : 



if t(peak) < 4min : p4 - 



t(peak) 



t(peak) 



where t(peak) is the time point of the peak of the pERK 
time response in minutes. 

5. The pERK time response curve should be smooth 
and not oscillate like an under-damped second order 
control system. If there were more than 3 peaks, the 
following penalty was defined: 



^5=n4 e «- 1 



i=l 



where m is equal to the number of FGF concentrations 
and d(i) is the number of peaks in pERK time response 
to FGF2 concentration L 



p2 



ERKfrac * totERK 
' max(pERK) 



where totERK is the number of molecules of total ERK 
present in the system, ERKfrac is the minimum fraction 



6. pERK dose-response curve should be biphasic at cer- 
tain time-points as observed in experiments. 



k / / .\ \ 2 

a{i) 
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where k is the number of time points, a(i) is the mimi- 
mal concentration of pERK at high FGF2 concentrations 
at time point i, ((i) is the approximate factor by which 
pERK should decrease from its maximal to its minimal 
response (determined from the experimental data) and 
is the maximal level of pERK for time point I C was 
chosen to be equal to [1.0, 1.0, 1.0, 1.0, 1.0, 0.95, 0.85, 
0.5, 0.5, 0.9, 0.9] for the different time points. 

The final objective function value was calculated as 
the product of all penalty values and the RMSE error. 

We compared the two objective functional forms for 
their ability to estimate the parameters of the FGFR 
model. The estimation was done 40 times, 20 with RMSE- 
based and 20 with hybrid objective function. For each of 
the 40 estimations, initialization was made completely 
randomly based on the range of parameter values pro- 
vided as input to the algorithm. The following key obser- 
vations were made - 

a. Although the final objective function value obtained 
was, in general, lower for the RMSE-based objective 
function than the hybrid objective function, param- 
eter estimates from RMSE-based function provided 
physically meaningless solutions. For instance, 19 
out of 20 parameter estimates predicted that less 
than 2% of total ERK gets phosphorylated by FGFR 
pathway. In contrast, 12 out of 20 parameter esti- 
mates obtained using hybrid objective function satis- 
fied all the experimental and physicality constraints. 
Each of these fits can be used as good initial esti- 
mates for obtaining more refined parameter 
estimates. 

b. Visual inspection of the fits showed that many of the 
parameter sets obtained using RMSE-based function 
showed only minimal to no biphasic dose response. 
Further investigation revealed that these parameter 
sets fit the low and intermediate doses of FGF2 ex- 
tremely well but completely failed to fit the response 
to high doses. In addition, many fits failed to match 
the time point of peak phosphorylation. 

Based on these observations, we hypothesize that the 
objective function manifold that uses a combination of 
RMSE and constraints gets rid of various local minima 
that correspond to physically meaningless parameter 
estimates and thus facilitates a much more efficient es- 
timation of parameter values. It can be expected that 
these findings can be extended to parameter estima- 
tion problems for various other signaling pathway 
models and other systems of ODEs. However, a formal 
comparison of the RMSE-based and hybrid objective 
function for benchmarked problems is beyond the 
scope of this work and remains to be verified in future 
studies. 



Section 4. Additional/modified equations and parameters 
for the heparin addition experiment (Table 4) 



dt 



ds2 
dt 

ds^ 
dt 



ds\Q 
dt 



ds\^ 
dt 



ds 23 
dt 

dS24 _ 

dt 

ds 25 
dt 



ds 2 7 
dt 

ds 28 
dt 

ds 29 
dt 

dsso 
dt 

ds 3 i 
dt 



"ft/la * S 2 * Si-krla * S 9 )-(kf 5a * Si * S 8 -k r5a * S W ) 
-ft/5«l * Si * s 24 -k r5a * S 25 ) 



- (kf 0 * s 2 * s 3 -k r0 * s 8 ) - (k/ia * s 2 * si-k 

- (kfQH * S 2 * S 23 -k r 0H * 524) 



-ft 



■r\ * S 9 ) 



C/15 * 5 4 * Sn-k r i5 * S 20 ) + kfdpl * 5i4 
~{kf43 * 519 * S4-k r 43 * Stf)" ft/15 * 5 4 * S 29 -/c r i 5 * S 32 ] 
~ ft/15 * H * 5 28 -/Crl5 * S33) 

= (kfsa * si * s 8 -k r5a * 5i 0 )-2 

* (kfdim * SiQ * SiQ-krdim * ^ll) 
-(kfdim * 5 25 * SiQ-krdim * ^26) 

= kf 19 * S 20 -kf dp l * 514- ft/35 * 5i4 * 5 5 -/r r35 * 5 2 l) 
+ kfse * 5 2 1- ft/43 * 5l9 * 5i 4 -/<r43 * 5l 8 ) 
+ kf 19 * 5 32 + kf 19 * 5 33 

= - ft/0// * 5 2 * S 23 -k r0 H * 5 24 ) 

(kfQH * s 2 * s 23 -k r0H * s 2 ^)-(kf 5al * si * s 2 4-k r5a * s 2 s) 

= ft/5«l * 5i * S 2 4-k r 5a * 5 2 5)~2 

* (kfdim * S 2 5 * S 2 5-k r dim * ^27) 
-(kfdim * 525 * SiQ-krdim * 526) 

= (kfdim * 525 * Sio-k r dim * S 26 )-I<jph * S 2 6 



1 (kfdim * 525 * S 25 -k r dim * S 27 )-I<jph * S 2 7 

: kfph * S 26 -(kfi n tl * S 2 s-k r i n tl * S 30 ) 

-ft/15 * 5 4 * S 28 -/C r i 5 * 533) + kf 19 * S33 

: kfph * S 27 -(kf int i * S 2 9-k rint i * S31) 

- ft/15 * S4 * s 29 -k r is * s 32 ) + kf 19 * s 32 
1 -(kfintl * S28-krintl * S30) 

= '(kfintl * s 29 -k rint i * S31) 
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Table 4 List of additional species and parameters in heparin perturbation model 



Name 


Species 














s23 


Heparin (HE) 














s24 


FGF:HE 














s25 


FGF:HE:FGFR 














s26 


(FGF:HE:FGFR) 2 














s27 


(FGF:HE:FGFR):(FGF:HSGAG:FGFR) 














s28 


n rfFGF-HF-FGFRV(FGF-HSGAG-FGFRY) 














s29 


|J\J \3\ .1 II 1 \J\ ly 2 














s30 


i((hGh:Hb:hGhK):(hGh:HbGAG:rGrK)j 














s31 


i(FGF:HE:FGFR) 2 














s32 


p(FGF:He:FGFR) 2 :FRS2 














s33 


p((FGF:HE:FGFR):(FGF:HSGAG:FGFR)):FRS2 














Name 


Unit Setl 


Set2 


Set3 


Set4 


Set5 


Set6 


Set7 


kfSal 


1/((molecule/cell)*s) 6.74E-1 1 


1.34E-10 


1.93E-10 


4.42E-1 1 


9.33E-1 1 


5.14E-10 


2.42E-08 


kfOH 


1/((molecule/cell)*s) 2.67E-08 


1.10E-08 


8.61 E-09 


1 .97E-08 


1 .46E-08 


1 .52E-08 


2.60E-08 


k r 0H 


1/s 1.94E-03 


8.00E-04 


6.29E-04 


1 .44E-03 


1 .07E-03 


1.11E-03 


1.95E-03 



= (*/15 * *4 * S 2 9-I<rl5 * S 32 ) ~kf\9 * 5 32 



ds 33 
dt 



[kfis * 5 4 * 5 2 8 -krlS * 533) ~kfl9 * 5 33 



Section 5: Additional equations and parameters for the 
MEKi addition experiment (Table 5) 



= ~{kf35 * S U * S 5 -k r35 * 5 2 l ) 

+ kfd p 2 * Si 5 -(kjMEKi * 5 23 * S 5 -k rM EKi * ^25) 



<i5 



^ = kf i 9 * S 2 0-kf dpl * 5i 4 -(/r /3 5 * 5i4 * 5 5 -/r r3 5 * 5 2 l ) 
+ kf 36 * S 2 l-(kf43 * 5i 9 * 5i 4 -/c r43 * 5i 8 ) 
- (kf 35 * $25 * 5i4-/<r35 * ^ + kf 36 * $26 



<i5 
dt 



J- = kf36 * S 2 \-kfdp2 * 5i 5 -(A:/39 * 5i 5 * 5 6 -/c r39 * 5 22 ) 
+ £/40 * S 2 2- (kjMEKi * ^23 * Sis-krMEKi * ^4) 



^5 23 



-{kfMEKi * ^23 * Sis-I< r MEKi * $24) 
-(kjMEKi * ^23 * Ss-krMEKi * ^5) 
-(kfMEKi * ^21 * S 2 3-k r MEKi * ^26) 

^24 /. 7 \ 

= [KjMEKi * ^23 * s \5~ KrMEKi * *24 J 

+ (AT/36 * S 2 6-kjdp3 * 524) 

= [KjMEKi * ^23 * Ss~ K r MEKi * $25 j 

+ (*/35 * 5 25 * 5i4-/<r35 * ^ + ^3 * *24 

ds ™ (1 1 \ 

= [KjMEKi * ^21 * S23- KrMEKi * S 2 6 ) 

+ (kf 35 * 5 2 5 * Su~I<r35 * S 2 6)~kf36 * *26 



Table 5 List of additional species and parameters in MEKi perturbation model 

Name Species 

s 23 MEKi (=U0126) 

s 24 pMEK:MEKi 
s 25 MEK:MEKi 
s 26 pFRS2:MEK:MEKi 



Name 


Unit 


Setl 


Set2 


Set3 


Set4 


Set5 


Set6 


Set7 


kfMEKi 


1/((molecule/cell)*s) 


2.61E-05 


4.14E-08 


3.28E-07 


8.1 5 E-06 


1 .00E-03 


8.73E-05 


5.78E-06 


krMEKi 


1/s 


1 .66E-06 


4.91 E-06 


6.60E-05 


3.94E-05 


2.99E-01 


3.55E-04 


1.04E-04 
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Section 6: Analyzing parameter uncertainty using 
Markov-chain Monte Carlo sampling. 

Parameter uncertainty for the seven sets identified using 
PSO was analyzed using Markov-chain Monte Carlo 
(MCMC) sampling approach as described previously 
[27]. Hug et. al. showed that a multi-chain sampling ap- 
proach can efficiently sample the posterior probability 
distribution for high-dimensional and non-linear models 
such as the FGF signaling model presented in this paper. 
Initializing from each of the seven parameter sets, 
Markov chains of length 100,000 were run using a modi- 
fied version of adaptive metropolis' algorithm as de- 
scribed in Haario et. al. [50]. Additional file 2: Figure S2 
shows the results of MCMC sampling for each of the 39 
parameters with the initial positions marked as red stars. 

For some parameters like 'H' or 'Kdla' , the estimated op- 
timal value for all seven optimal parameter sets are close 
together and located at the boundary of the allowed par- 
ameter range. Accordingly, the MCMC samples show par- 
ameter distribution close to the boundary. Importantly, 
the distribution shows that other parameter sets with 
values away from the boundary are able to explain the ex- 
perimental data equally well. Thus, this analysis gives a 
more realistic impression of the objective function mani- 
fold as a function of parameter values. For some other 
parameters like l<f44', two distinct clusters of MCMC sam- 
ples were observed. This indicates that the parameter sets 
belong to two distinct basins of attraction. In summary, 
these results indicate that the seven sets of parameters 
identified by PSO are representative for the distribution 
of good fits that can be identified by MCMC sampling 
and provide the argument for limiting the further model 
analysis to these representative sets. The model was fur- 
ther analyzed using this large set of parameter values. 
For each of the parameter set represented in Additional 
file 2: Figure S2, the model satisfies all the qualitative con- 
straints as described in Materials and methods section 
3, including the constraint for biphasic pERK response. 
Thus, the model predicts biphasic pERK response for a 
wide range of parameter values represented in Additional 
file 2: Figure S2. 



Additional files 



Additional file 3: Figure S3. Detailed schematic of the model-reaction 
network for FGFR pathway. 
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